Xi-Language Reference: Library Functions

Library functions are a bunch of functions written in Xi that sometimes come handy. They are also good demos of the capabilities of Xi.

fft_image_filter (filters the spectrum of an image)

/***********************************************************************
*Name: fft_image_filter
*Type: BuildinFuntion
*Section: 62 Geometrical Transformations
*Parameters: char[] image, int radius
*Return: char[] (contains the filtered image)
*DescriptionShort: filtering of an image in the frequency domain
*DescriptionLong:
 The function \e\fft_image_filter\e\ performs an ideal low pass
 filter with cutoff frequency \e\radius\e\. 
*Example:
*SeeAlso: fft, rev_fft
*End:
************************************************************************/

char[] fft_image_filter(char image[], int radius)
{
   s=size(image);
   rows=s[3];
   cols=s[4];

   mrows=rows/2;
   mcols=cols/2;

   b=carr(rows,cols);
   b[mrows-radius:mrows+radius, mcols-radius:mcols+radius]=1;
   b=shift(b,mrows,mcols);

   n=abs(rev_fft(fft(image)*b));

   delete b;
   delete s;

   return (char) (n / (max(n)/255));
}


interval (Subdivides an interval)

/***********************************************************************
*Name: interval
*Type: BuildinFunction
*Section: 515 Array functions
*Parameters: (double) x0, (double) xn, (int) n = 100
*Return: double[] (contains the x-values)
*DescriptionShort: Subdivid a interval into equidistant tiles
*DescriptionLong:
 \t\interval\t\ simply subdivides the interval [\e\x0\e\,\e\xn\e\] into
 \e\n\e\ equistant tiles. It returns an array with \e\n\e\+1 elements, such
 that the first element is \e\x0\e\ and the last is \e\xn\e\.
*Example:
 (  1)>x=interval(0,1,10);
 Function interval defined
 (  2)>print(x)
 
            0          0.1          0.2          0.3          0.4          0.5
          0.6          0.7          0.8          0.9            1
*SeeAlso: interval2d
*End:
************************************************************************/
 
double [] interval(double x0, double xn, int n=100)
{
  return dincarr(n+1)/n*(xn-x0)+x0;
}

interval2d (Subdivides a rectanglular region)

***********************************************************************
*Name: interval2d
*Type: BuildinFunction
*Section: 515 Array functions
*Parameters: (double) x0, (double) y0, (double) xn, (double) yn,
             (int) nx = 100, (int) ny = 100
*Return: [double[],double[]] (contains the x- and y-values)
*DescriptionShort: Subdivid a rectangle into tiles
*DescriptionLong:
 \t\interval2d\t\ generates a tiling of a rectangular domain of the plane.
 This is especially useful for evaluation of two dimensional functions.
 The rectangle is defined by the two corners (\e\x0\e\,\e\y0\e\) and
 (\e\xn\e\,\e\yn\e\). \t\interval2d\t\ will return two
 (\e\nx\e\+1)x(\e\ny\e\+1) array contains the x- resp. y-coordinates,
 such that x[.,0]=\e\x0\e\, x[.,\e\nx\e\]=\e\xn\e\ , y[0,.]=\e\y0\e\ and
 y[\e\ny\e\,.]=\e\yn\e\.
*Example:
 (  1)>[x,y]=interval2d(0,0,1,1,2,2);
 Function interval2d defined
 (  2)>print(x,y)
 
            0          0.5            1
            0          0.5            1
            0          0.5            1
 
            0            0            0
          0.5          0.5          0.5
            1            1            1
*SeeAlso: interval2d
*End:
************************************************************************/
  
[ double[], double[] ] interval2d(double x0, double y0,
                                  double xn, double yn,
                                  int nx=100, int ny=100)
{
  x=y=darr(nx+1,ny+1);
  for(i=0;i<=nx;i++)
    for(j=0;j<=ny;j++)
      {
        x[i,j]=i*(xn-x0)/nx+x0;
        y[i,j]=j*(yn-y0)/ny+y0;
      }
  return [x,y];
}

rainbow(generates a rainbow colormap)

/************************************************************************
*Name: rainbow
*Type: Xi function
*Section: 61 Color Manipulations
*Parameters: (int) n = 255, (double) s = 0.8, (double) v=0.8
*Return: char[] (contains the colormap)
*DescriptionShort: generates a rainbow colormap
*DescriptionLong:
The function \e\rainbow\e\ produces a rainbow colormap with \e\n\e\ entries,
the saturation (determines how pure a color is) \e\s\e\ and the
value (intensity of the color) \e\v\e\.
*Example:
>rainbowMap=rainbow();
Function rainbow defined
>loadct(rainbow(126,0.4,0.4));
>window(0,\bpp=7);
>show(cincarr(126,1)+2);
*SeeAlso: graymap, loadct
*Reference:
*End:
*************************************************************************/

char[] rainbow(int n=255, double s=0.8, double v=0.8)
{
   double hsv[n,3];
   phi_step=280.0/n;

   for (i=0; i


velocity_field (draw a velocity field)

int arrowHead(double x,  double y,
              double nx, double ny,
              double headLen)
{
   dx=nx-x;
   dy=ny-y;

   vc={dx, dy};
   r=sqrt(dx*dx+dy*dy);
   if (r == 0)
     sc=0.0;
   else 
     sc=headLen/r;

   D={ { 0.866025 , 0.5 }, { -0.5, 0.866025 } };
   dn1=D#vc           *sc;
   dn2=transpose(D)#vc*sc;

   line(nx,ny,nx-dn1[0],ny-dn1[1]);
   line(nx,ny,nx-dn2[0],ny-dn2[1]);

   return 0;
}

/************************************************************************
*Name: velocity_field
*Type: Xi function
*Section: 521 2D Plots
*Parameters: (double[]) u, (double[]) v,
             (double[]) x = -1, (double[]) y = -1,
             (int) nstep = 1, (double) stepSize = 0.1, 
             (double) headScale = 0.5, 
             (int) noaxis = 0, (int) curve = 0
*Return: -1
*DescriptionShort: draws a velocity_field
*DescriptionLong: \e\velocity_field\e\  plots a velocity field with arrows 
 pointing in the
 direction of the given velocity field. \e\u\e\ and \e\v\e\ give the
 X and Y components of the velocity field.

 If the parameters \e\x\e\ and \e\y\e\ are given the should be arrays.
 The size of \e\y\e\ must be equal to the first dimension of u, v and x
 must be equal to the secound dimension of u, v. 
 \e\startX\e\ and \e\startY\e\ set the start points and the
 array \e\color\e\ sets the colors  of each arrow of each arrow.
 \e\nsteps\e\ changes the number of segments of each arrow and
 \e\stepSize\e\ determine the number of segments.
 The parameter \e\headScale\e\ is the proportional factor between
 the length of the arrow and sides of the arrow-head.
 Corresponding to the \e\plot\e\ command the flag \e\noaixs\e\
 suppresses the axis.
 If the segments of the arrows should be connected with splines the
 flag \e\curve\e\ will do the job for you.
*Example:
 >u=cos(dincarr(10,15)/20.);
 >v=transpose(sin(dincarr(15,10)/20.));
 >contour(q,\nlevel=7,\curve);
*SeeAlso: velocity_field3d
*Reference:
*End:
*************************************************************************/


int velocity_field(double u[], double v[],
                   double x[]     ={-1}, double y[]     ={-1},
                   double startX[]={-1}, double startY[]={-1},
                   int    nsteps   =1,
                   double stepSize =0.1,
                   double headScale=0.5,
                   int noaxis=0, int curve=0)
{
   nsteps++;
   s1=size(u);
   s2=size(v);
   x1=s1[4]; x2=s2[4];
   y1=s1[3]; y2=s2[3]; 

   if (s1[1] != 2 || s2[1] != 2)
    {
      error("velocity_field","Arrays u, v must have two dimensions");
      return -1;
    }
   if (x1 != x2 || y1 != y2)
    {
      error("velocity_field","Arrays u, v  must have the same range");
      return -1;
    }

   s3=size(x);
   if (s3[2+s3[1]] == 1)
     {
       x=dincarr(x1);
       s3=size(x);
     }
   else
     if (s3[1] != 1 || s3[3] != x1)
       {
         error("velocity_field","Array x have a wrong dimension or range");
         return -1;
       }

   s4=size(y);
   if (s4[2+s4[1]] == 1)
     {
       y=dincarr(y1);
       s4=size(y);
     }
   else
     if (s4[1] != 1 || s4[3] != y1)
       {
         error("velocity_field","Array y have a wrong dimension or range");
         return -1;
       }

   s5=size(startX);
   if (s5[2+s5[1]] == 1)
     {
       startX=darr(x1*y1);
       for (i=0; i= x[x1-1])
               { aj=x1-2; t=1.0; }

           for (m=0; m < y1-1; m++)
             if (y[m] <= py[k-1] && py[k-1] < y[m+1])
               {
                 ak=m;
                 s=(py[k-1]-y[ak])/(y[ak+1]-y[ak]);
               }

	   if (py[k-1] < y[0])
             { ak=0; s=0; }
           else
             if (py[k-1] >= y[y1-1])
               { ak=y1-2; s=1; }

           dpx[k]= (1-s)*(1-t)*dx[ak,aj] + s*(1-t)*dx[ak+1,aj]
                  +s*t*dx[ak+1,aj+1]     + (1-s)*t*dx[ak, aj+1];
           dpy[k]= (1-s)*(1-t)*dy[ak,aj] + s*(1-t)*dy[ak+1,aj]
                  +s*t*dy[ak+1,aj+1]     + (1-s)*t*dy[ak, aj+1];

           px[k]=px[k-1]+dpx[k];
           py[k]=py[k-1]+dpy[k];	    
          }

       if (curve)
         plot(px,py,\curve);
       else
         plot(px,py,\line);

       arrowHead(px[nsteps-2], py[nsteps-2],
                 px[nsteps-1], py[nsteps-1],
                 headScale*total(sqrt(dpx*dpx+dpy*dpy))); 
     }

   endbuffer();
   return 0;
}

velocity_field3d (draws a 3 dimensional velocity field)

int arrowHead3d(double x,  double y, double z,
                double nx, double ny, double nz,
                double headLen,
                char color=1)
{

   dx=nx-x;
   dy=ny-y;
   dz=nz-z;

   vc={dx, dy};
   r=sqrt(dx*dx+dy*dy);
   if (r == 0)
     sc=0.0;
   else 
     sc=headLen/r;

   D={ { 0.866025 , 0.5 }, { -0.5, 0.866025 } };
   dn1=D#vc           *sc;
   dn2=transpose(D)#vc*sc;
   dn3=dz*sc;

   line3d(nx,ny,nz,nx-dn1[0],ny-dn1[1],nz-dn3,\color=color);
   line3d(nx,ny,nz,nx-dn2[0],ny-dn2[1],nz-dn3,\color=color);

   return 0;
}


/************************************************************************
*Name: velocity_field3d
*Type: Xi function
*Section: 522 3D Plots
*Parameters: (double[]) u, (double[]) v, (double[]) w,
             (double[]) x = -1, (double[]) y = -1, (double[]) z = -1,
             (double[]) startX = -1,
             (double[]) startY = -1, (double[]) startZ = -1,
             (char) color[] = -1 
             (int) nstep = 1, (double) stepSize = 0.1, 
             (double) headScale = 0.5, 
             (int) noaxis = 0
*Return: -1
*DescriptionShort: draws a velocity_field
*DescriptionLong: \e\velocity_field3d\e\  plots a velocity field with 
 arrows pointing in the direction of the given velocity field. 
 \e\u\e\, \e\v\e\ and \e\w\e\ give the X, Y and Z components of the
 3 dimensional velocity field.
 If the parameters \e\x\e\, \e\y\e\ and \e\z\e\ are given the should
 be arrays. The size of \e\y\e\ must be equal to the first
 dimension of u, v, w and x must be equal to the second dimension
 of u, v, w. At least the size of \z\  
 must be equal to the third dimension of u, v, w.
 \e\startX\e\, \e\startY\e\ and \e\startZ\e\ set the start points
 and the array \e\color\e\ sets the colors  of each arrow.
 \e\nsteps\e\ changes the number of segments of each arrow and
 \e\stepSize\e\ determine the number of segments.
 The parameter \e\headScale\e\ is the proportional factor between
 the length of the arrow and sides of the arrow-head.
 Corresponding to the \e\plot\e\ command the flag \e\noaixs\e\
 suppresses the axis.
*Example:
*SeeAlso: velocity_field
*Reference:
*End:
*************************************************************************/


int velocity_field3d(double u[], double v[], double w[],
                     double x[]     ={-1},
                     double y[]     ={-1},
                     double z[]     ={-1},
		     double startX[]={-1},
                     double startY[]={-1},
                     double startZ[]={-1},
                     char   color[] ={-1},
                     int    nsteps   =1,
                     double stepSize =0.1,
                     double headScale=0.5,
                     int noaxis=0)
{
   nsteps++;
   su=size(u);
   sv=size(v);
   sw=size(w);

   xu=su[5]; xv=sv[5]; xw=sw[5];
   yu=su[4]; yv=sv[4]; yw=sw[4];   
   zu=su[3]; zv=su[3]; zw=sw[3];

   xl=su[5];
   yl=su[4];
   zl=su[3];

   if (su[1] != 3 || sv[1] != 3 || sw[1] != 3)
    {
      error("velocity_field","Arrays u, v and w must have two dimensions");
      return -1;
    }
   if (xu != xv || xu != xw || yu != yv || yu != yw || zu != zv || zu != zw )
    {
      error("velocity_field","Arrays u, v and w  must have the same range");
      return -1;
    }

   s3=size(x);
   if (s3[2+s3[1]] == 1)
     {
       x=dincarr(xl);
       s3=size(x);
     }
   else
     if (s3[1] != 1 || s3[3] != xl)
       {
         error("velocity_field","Array x have a wrong dimension or range");
         return -1;
       }

   s4=size(y);
   if (s4[2+s4[1]] == 1)
     {
       y=dincarr(yl);
       s4=size(y);
     }
   else
     if (s4[1] != 1 || s4[3] != yl)
       {
         error("velocity_field","Array y have a wrong dimension or range");
         return -1;
       }

   s5=size(z);
   if (s5[2+s5[1]] == 1)
     {
       z=dincarr(zl);
       s5=size(z);
     }
   else
     if (s5[1] != 1 || s5[3] != yl)
       {
         error("velocity_field","Array z have a wrong dimension or range");
         return -1;
       }

   s6=size(startX);
   s7=size(startY);
   s8=size(startZ);

   if (s6[2+s6[1]] == 1)
     {
       startX=darr(xl*yl*zl);
       for (i=0; i= x[xl-1])
               { aj=xl-2; t=1.0; }

	   for (m=0; m < yl-1; m++)
             if (y[m] <= py[k-1] && py[k-1] < y[m+1])
               {
                 ak=m;
                 s=(py[k-1]-y[ak])/(y[ak+1]-y[ak]);
               }

	   if (py[k-1] < y[0])
             { ak=0; s=0; }
           else
             if (py[k-1] >= y[yl-1])
               { ak=yl-2; s=1; }


          for (m=0; m < zl-1; m++)
             if (z[m] <= pz[k-1] && pz[k-1] < z[m+1])
               {
                 al=m;
                 r=(pz[k-1]-z[al])/(z[al+1]-z[al]);
               }

	   if (pz[k-1] < z[0])
             { al=0; s=0; }
           else
             if (pz[k-1] >= z[zl-1])
               { al=zl-2; s=1; }


           dpx[k]=  (1-s)*(1-t)*(1-r)*dx[al,ak,aj]
                  + s*(1-t)*(1-r)    *dx[al+1, ak, aj]
                  + (1-s)*t*(1-r)    *dx[al, ak+1, aj]
                  + (1-s)*(1-t)*r    *dx[al, ak, aj]
                  + (1-s)*t*r        *dx[al, ak+1, aj+1]
                  + s*(1-t)*r        *dx[al+1, ak, aj+1]
                  + s*t*(1-r)        *dx[al+1, ak+1, aj]
                  + s*t*r            *dx[al+1, ak+1, aj+1];

           dpy[k]=  (1-s)*(1-t)*(1-r)*dy[al,ak,aj]
                  + s*(1-t)*(1-r)    *dy[al+1, ak, aj]
                  + (1-s)*t*(1-r)    *dy[al, ak+1, aj]
                  + (1-s)*(1-t)*r    *dy[al, ak, aj]
                  + (1-s)*t*r        *dy[al, ak+1, aj+1]
                  + s*(1-t)*r        *dy[al+1, ak, aj+1]
                  + s*t*(1-r)        *dy[al+1, ak+1, aj]
                  + s*t*r            *dy[al+1, ak+1, aj+1];

           dpz[k]=  (1-s)*(1-t)*(1-r)*dz[al,ak,aj]
                  + s*(1-t)*(1-r)    *dz[al+1, ak, aj]
                  + (1-s)*t*(1-r)    *dz[al, ak+1, aj]
                  + (1-s)*(1-t)*r    *dz[al, ak, aj]
                  + (1-s)*t*r        *dz[al, ak+1, aj+1]
                  + s*(1-t)*r        *dz[al+1, ak, aj+1]
                  + s*t*(1-r)        *dz[al+1, ak+1, aj]
                  + s*t*r            *dz[al+1, ak+1, aj+1];

           px[k]=px[k-1]+dpx[k];
           py[k]=py[k-1]+dpy[k];
           pz[k]=pz[k-1]+dpz[k];
          }

         plot3d(px,py,pz,\line,\color=color[i]);
	 arrowHead3d(px[nsteps-2], py[nsteps-2], pz[nsteps-2],
                     px[nsteps-1], py[nsteps-1], pz[nsteps-1],
                     headScale*total(sqrt(dpx*dpx+dpy*dpy)),
                     color[i]);
      }

   endbuffer();
   return -1;
}